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ABSTRACT 

We present a new non-parametric method for determining mean 3D density 
and mass profiles from weak lensing measurements around stacked samples of 
galaxies or clusters, that is, from measurement of the galaxy-shear or cluster-shear 
correlation functions. Since the correlation function is statistically isotropic, 
this method evades problems, such as projection of large-scale structure along 
the line of sight or halo asphericity, that complicate attempts to infer masses 
from weak lensing measurements of individual objects. We demonstrate the 
utility of this method in measuring halo profiles, galaxy-mass and cluster-mass 
cross-correlation functions, and cluster virial masses. We test this method on 
an N-body simulation and show that it correctly and accurately recovers the 3D 
density and mass profiles of halos. We find no evidence of problems due to a mass 
sheet degeneracy in the simulation. Cross-correlation lensing provides a powerful 
method for calibrating the mass-observable relation for use in measurement of 
the cluster mass function in large surveys. It can also be used on large scales to 
measure and remove the halo bias and thereby provide a direct measurement of 
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1. Introduction 

Gravitational lensing has become a powerful tool for studying cosmology and especially 
for studying the unseen matter which dominates the Universe. Weak gravitational lensing 
produces an image of the Universe which is slightly distorted — these distortions reveal the 
quantity and clustering properties of the elusive dark matter. One major application of weak 
lensing has been measurement of the mass distribution around rich galaxy clusters (Fahlman 
et al. 1994; Tyson & Fischer 1995; Luppino & Kaiser 1997; Clowe et al. 1998; Joffre et al. 
2000; Irgens et al. 2002; Cypriano et al. 2004) This has allowed the reconstruction of 2D 
projected cluster mass density maps which, combined with optical and X-ray measurements, 
has enabled measurement of mass-to-light ratios and determination of the cluster baryon 
fraction. 

Another important use of weak lensing is the measurement of the power spectrum of 
mass fluctuations through cosmic shear, which has now been detected by several groups (Van 
Waerbeke et al. 2000; Wittman et al. 2000; Bacon et al. 2000; Maoli et al. 2001; Rhodes 
et al. 2001; Hoekstra et al. 2002; Jarvis et al. 2003; Bacon et al. 2003; Hamana et al. 2003). 
Cosmic shear, or shear-shear correlations, probes cosmology through its dependence on both 
the growth of structure and the expansion history of the Universe (Hu 2002; Abazajian & 
Dodelson 2003; Takada & Jain 2004; Song & Knox 2004; Refregier et al. 2004) since it is 
sensitive only to the matter distribution, cosmic shear is theoretically clean. 

While cluster lensing allows the study of the mass distribution of individual objects, 
cosmic shear is a purely statistical quantity analogous to the angular power spectrum of 
the cosmic microwave background. There is another method, called galaxy-galaxy lensing 
or galaxy-shear correlations, which is a hybrid of the two. Galaxy-galaxy lensing measures 
the mean weak lensing shear produced by a sample of lens galaxies by averaging the shear 
signal over large numbers of them. It is thus a statistical measure which can be applied to 
different samples of galaxies to measure the average mass distribution around different galaxy 
types. After early null results (Kristian 1967; Tyson et al. 1984), the first detections using 
this technique were made by Brainerd et al. (1996); dell'Antonio & Tyson (1996); Griffiths 
et al. (1996); these measurements have quickly improved (Hudson et al. 1998; Fischer et al. 
2000; Wilson et al. 2001; McKay et al. 2001; Smith et al. 2001; Hoekstra et al. 2003). The 
Sloan Digital Sky Survey (SDSS) (York et al. 2000) has provided an excellent sample for 
measuring galaxy-galaxy lensing: Sheldon et al. (2004) measured the weak lensing signal 
around the SDSS main spectroscopic sample to high accuracy (see also (Seljak et al. 2005)). 
The galaxy-galaxy lensing signal is best interpreted, in the language of large-scale structure, 
as an estimate of the galaxy-mass correlation function. Combined with the galaxy auto- 
correlation function, it gives us information about galaxy biasing; on large scales, it provides 
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a direct measurement of <7gQ m that assumes only the physics of general relativity. 

Another arena in which such statistical weak lensing methods can be applied is the 
study of galaxy clusters. The number counts of galaxy clusters as a function of mass and 
redshift is a sensitive probe of cosmological parameters. The cluster mass function has been 
measured by several groups using various proxies for cluster mass (Bahcall & Cen 1992; 
White et al. 1993; Viana & Liddle 1999; Reiprich & Bdhringer 2002; Bahcall et al. 2003). 
The measurement of cluster number counts at low redshift is usually used to constrain the 
combination cr 8 f2^ with (3 ~ 0.5, although the exact value of j3 really depends on which range 
of mass scales is being probed (Rozo et al. 2004). Measurement of the evolution of the mass 
function with redshift allows one to break the degeneracy between these two parameters 
(Bahcall & Bode 2003). Galaxy clusters are strongly biased compared to the dark matter, 
and this bias as a function of halo mass is a prediction of theoretical models (Kaiser 1984; 
Mo et al. 1996; Sheth & Tormen 1999; Seljak & Warren 2004). A measurement of cluster 
bias as a function of mass via cluster correlations thus provides an independent constraint 
on erg (Rozo et al. 2004). In the future, planned large-area surveys will aim to determine the 
properties of the dark energy (its equation of state and energy density) by counting clusters 
and measuring their clustering to redshifts z ~ 1 — 2 (Haiman et al. 2001). 

One of the main challenges to using galaxy clusters to constrain cosmology is uncertainty 
in the determination of clusters masses. Typically one selects clusters and/or infers their 
masses based on some observable proxy for mass, such as X-ray temperature or luminosity 
(Borgani & Guzzo 2001; Del Popolo 2003), Sunyaev-Zel'dovich effect (SZE) flux decrement 
(Grego et al. 2001), optical richness, or galaxy velocity dispersion (van der Marel et al. 2000; 
McKay et al. 2002; Lokas & Mamon 2003). In all these cases, the mass-observable relation, 
as well as its scatter and possible evolution with redshift, must be calibrated (Pierpaoli et al. 
2003) either by simulations, by combining number counts with spatial clustering information 
(Majumdar & Mohr 2004; Lima & Hu 2004, 2005), or by 'direct' mass measurements of a 
subsample of the clusters. X-ray and SZE measurements are sensitive to varying degrees to 
the physics of the intra-cluster medium, such as radiative transfer, cooling, accretion, and 
feedback mechanisms associated with star and galaxy formation, possibly leading to scatter 
in the mass-observable relation. Dynamical mass measurements through the galaxy velocity 
dispersion suffer from uncertainties in the anisotropy of the velocity distribution function, 
possible bias of the galaxy velocities relative to the dark matter, and possible violation of 
the required assumption of dynamical equilibrium (e.g., due to recent mergers). 

Weak lensing provides a method for measuring cluster masses directly and thereby 
calibrating mass-observable relations. The advantage of lensing is that it is sensitive to the 
total mass regardless of its makeup or dynamical state; the only physics involved is the 
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deflection law given by general relativity. However, lensing measurements have limitations 
of their own. Weak lensing mass reconstruction of individual clusters can only obtain a 
two-dimensional surface mass density; information about the structure along the line of 
sight is irretrievably lost. Typically these measurements come from pointed observations 
of known massive clusters, selected by X-ray luminosity or other property. Because the 
observations need to go deep enough to obtain the necessary density of background source 
galaxies to make spatially resolved mass maps, they are time-consuming and therefore tend 
to be limited in area, usually not extending much beyond the virial radius. Because of this 
limited range of physical scales covered, there are uncertainties arising from the mass sheet 
degeneracy (Bradac et al. 2004). In addition, large-scale structure along the line of sight 
and chance projections with clusters at different redshifts can influence individual cluster 
mass determinations, causing ~ 30% scatter around the true mass and possibly a bias in 
the mass determination (Cen 1997; Metzler et al. 1999; Hoekstra 2001; Metzler et al. 2001; 
White et al. 2002; Hoekstra 2003; Clowe et al. 2004; Dodelson 2004). 

In this paper, we demonstrate a new way to measure the average density and mass 
profile of clusters using stacked samples that, in principle, can be selected by any observable 
proxy. We introduce a completely non-parametric method to deproject the average shear 
profile and obtain the mean 3D density and mass profiles for a sample. This method only 
assumes the statistical isotropy of the cluster-mass correlation function. We perform tests 
on an N-body simulation to show that the method correctly recovers the average 3D density 
and mass profiles. From the derived mass profile of the clusters, one can measure a virial 
radius and virial mass without assuming a model for the profile. Unlike individual cluster 
lensing reconstructions, these mass estimates are not affected by large-scale structure or 
chance projections along the line of sight. We show that these methods will allow one to 
constrain the density profiles of galaxy clusters, calibrate mass-observable relations, and 
thereby measure the cluster mass function and bias. 

In §2, we derive the relations between tangential shear measurements and the 2D surface 
mass density. In §3, we invert the relation to obtain the 3D density profile. In §4, we present 
methods of deriving 3D mass profiles, and in §5 we describe some of the practical steps 
(interpolation, extrapolation, and error propagation) needed to implement the method. In 
§6, we apply the method to a large N-body simulation, showing that it accurately recovers 
the correct 3D density and mass profiles statistically, and in §7 we discuss why the method is 
insensitive to asphericity of halos and projected large-scale structure along the line of sight. 
§8 presents our conclusions and outlook for the future. In future papers, we will apply this 
method to measurement of the cluster-mass correlation function for optically selected clusters 
in the SDSS and in N-body simulations, derive scaling relations between optical richness and 
virial mass, and constrain the bias and the mass power spectrum by comparison with the 
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cluster-cluster correlation function. 



2. 2D Density Inversion 

The primary measurement in cross-correlation lensing is of the average tangential dis- 
tortion of background source galaxy images by foreground lens objects, typically galaxies or 
clusters. This distortion is measured in several radial annuli around the center of the lens 
and repeated for all lenses in the sample. For the purpose of this paper, we assume that 
we have redshift measurements for the lenses, as is the case with the SDSS and some other 
current and planned surveys. We will also assume that we have an unbiased estimate of the 
lensing strength, < Sj r * t >, which is to say we have a good understanding of the redshift 
distribution of the source galaxies, possibly from using photometric redshifts. 

For any set of lenses, we can define the lens-mass correlation function, £ Jm , and write 
the average 3D mass density of the lens population as 

p(r)=p [l+& m (r)]=p + Ap(r), (1) 

where the cosmic mass density p = Q m p crit and Ap(r) = p £/ m (r). The lensing deflection 
and the shear are determined by the excess mass density projected along the line of sight, 
S(i?) = f dzAp(r), where we have dropped the constant term for reasons described below. 
Here and below, R denotes the 2D separation from the lens center in the lens plane, r denotes 
the 3D radial separation from the lens, and z is the distance from the lens along the line of 
sight, so that r 2 = R 2 + z 2 . The average tangential shear, 7 T , in a thin annulus of radius R 
is given by 

S c „ t lT (R) = £(< R) - = AE(R), (2) 
where E(< R) is the mean surface density interior to radius R, the critical surface density is 

c 2 D s 

and D s , Di, and Di s are the angular diameter distances to the source, to the lens, and from 
lens to source. Defining the convergence k(R) = E(i2)/E cr j t , the expression for the shear 
can be put in the dimensionless form 

<y T (R) = R(<R)-K(R). (4) 

Equation 2 defines AE(i?), the difference between the mean surface density inside a disk 
of radius R and E(i?) (azimuthally averaged), which is the primary observable in cross- 
correlation lensing. From this definition, one can see that a uniform-density mass sheet 
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of transverse size > R produces no shear; this is referred to as the mass sheet degeneracy 
(Bradac et al. 2004) and is the reason why we dropped the constant term in the definition 
of S(i?) above. 

We have shown how to predict the lensing observable AS from the theoretical quantity 
of interest, the correlation function Ap = £l m £i m . We would like to invert this relation, to 
infer fl m C,i m from measurements of AS. One way to proceed would be to use a parametric 
model for £j m , such as a power-law or NFW profile (Navarro et al. 1997). If a model profile 
is assumed, one can project the model to 2D, compute AS(i?), fit the shear data, and 
determine the best-fit model parameters. This method has the advantages of simplicity and 
of small formal error bars on the derived model parameters. Of course, the problem with any 
model-dependent method is that an incorrect model will result in incorrect interpretation 
of the best-fit parameters. In addition, parametric methods force smoothness on the data 
and thus do not allow one to see any features in £j m or detect fine-scale deviations from 
the model. Since we would like to constrain the full density profile and not just determine 
its best-fit parameters for a particular model, in this paper we develop a non-parametric 
approach which inverts the process above and allows a direct estimate of the 3D density 
profile Ap(R). We start by deriving a 2D reconstruction for the mean convergence profile, 



where ' denotes d/dR. This is just the azimuthally symmetric version of Kaiser's equation 
(Kaiser 1995) 



where 7, denotes the two-component shear field, and 7^ = d^i/ddj denotes the derivative 
in the lens plane. Equation 6 is the starting point for most finite-field non-linear inversion 
methods for weak lensing (Kaiser 1995; Lombardi & Bertin 1999; Seitz & Schneider 2001). 
Note that these equations relate the shear to the density in a local way, in contrast to 
Equation 4 for the tangential shear, which is non-local. As we will show below, both the 
local and non-local aspects of these equations are useful when inferring the mass profile. 

Schneider & Seitz (1995) pointed out that, except in the limit of very weak fields, 
7, k <C 1, the shear, 7, is not directly observable. The average shape of background objects 
is not determined by 7 or k but rather by the degenerate combination g = 7/(1 — k), 
which is sometimes referred to as the reduced shear. Furthermore, the average ellipticities 
of background galaxies is equal to the complex distortion, 5 = 2g/(l + \g\ 2 )- Solving this 



k(R). 



Taking the derivative of Equation 4 with respect to R one can show that 



k\R) = 1 ' t {R) + 2 1t {R)/R ) 



(5) 




(6) 
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quadratic equation for g results in g — (1 ± a/1 — \5\ 2 )/5*. Since there are two roots to this 
equation, one cannot tell from purely local measurements which root to choose. This reflects 
an inherent local degeneracy to any weak lensing measurements. For weak fields, one should 
take the negative root, for which 7 ~ g ~ 5/2. The choice of sign is related to image parity; 
it changes when crossing critical curves in the lens plane, where the magnification diverges 
and arcs are present. In practice, it is only in the cores of massive clusters where one might 
make the wrong choice of parity, and this is unlikely to present a major problem for the 
statistical method of weak lensing we pursue here. These rare, extreme regions are not the 
places where weak lensing will be the most useful probe of the mass distribution. 

The other continuous degeneracy is not as problematic: one can show that if g is mea- 
sured, one can invert for k up to the mass sheet degeneracy (Kaiser 1995). For our az- 
imuthally symmetric case, we can simply substitute g(l — k) for 7 in Equation 4, resulting 
in 

K ' d wi ^ g' + 2g/R _ n( r>\ C7\ 

= dR ln(1 - K) = ~T^g- = (7) 
so that ln(l — k) is determined up to a constant. Imposing the boundary condition k{oo) = 0, 
one can integrate this ordinary differential equation: 

PCX) 

k(R) = 1 - exp(- / dp G(p)), (8) 
Jr 

which in the weak-field limit becomes 

POO 

k(R)= / dp [ 7 t(p)' + 2 7 t(p)/p]. (9) 
Jr 

In section 5, we will discuss the truncation of integrals of this kind to a finite region. Since 
Equation 9 is linear in j T , one could simply multiply it by X crit to obtain an expression for 
£(.R). This is suitable for an analysis such as that in Sheldon et al. (2004), where AE = S crit 7 
is averaged. Beyond the weak field regime, the situation is more complicated, because of 
the 1 — g term in the denominator of G(R); in this case, an average of the lensing distortion 
does not simply yield an average of S(i?). Nonetheless, as noted above, this should only 
be a problem for small scales and high densities, i.e., for the cores of very massive clusters, 
and in that regime one could apply a correction or an iterative process to determine S(i?). 
Moreover, there is only a small spatial range where this correction is important and where 
one can still safely assume that 7 = g~5/2is the correct root. The method of averaging 
shear profiles may simply not be the most prudent method for probing the non-weak regime; 
strong lensing can likely tell us more about the density profiles of massive clusters at such 
small scales. Alternatively, other methods, such as the maximum likelihood approach of 
Schneider & Rix (1997), might be able to better handle this regime, though at the cost of 
introducing explicit model dependence. 
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3. Inverting to 3D Density Profiles 



We can rewrite Equation 7 (after multiplying by S cr jf) as 



E'(R) = E crit (1 - k)G(R) = Y> cnt G(R) exp(- / dp G(p)) 



(10) 



or in terms of the shear, (multiplying Equation 5 by Tl crit ) 



E'(R) = AE'(R) + ^AZ(R) 



(11) 



An equation for T,'(R) is useful because of the existence of the Von Zeipel (or sometimes 
called Abel) inversion formula for the 3D density profile: 



where, as before, Ap(r) = p(r) — p, because we do not recover the average mass density of 
the Universe: lensing is not sensitive to mass sheets. That is, we have arbitrarily chosen 
the boundary condition p(oo) = p or Ap(oo) = 0; for further discussion of the mass sheet 
degeneracy, see §5.2. 

Inversions of this type have a long history in astronomy. Von Zeipel (1908) derived 
this inversion formula and used it to determine the 3D density profile of globular clusters 
from imaging data. Von Zeipel's proof of this inversion formula rests on reducing it through 
substitutions to Abel's formula (Abel 1826), which in turn, is usually proven through the use 
of the Laplace transform. Plummer (1911), also concerned with globular clusters, derived 
a similar formula which uses intensities in long parallel strips rather than circular annuli. 
Recently Kaastra (1989) has shown that Von Zeipel's formula can be derived from Plum- 
mer's formula much more easily without using Abel's formula. The interested reader should 
consult Kaastra (1989) as well as Bremer (1995) and references therein. These inversions 
have been used for astronomical deprojections in many different contexts, including globu- 
lar clusters, galaxy clusters, elliptical galaxies, supernova remnants, planetary nebulae, and 
galaxy correlation functions. 




(12) 



The inverse of this equation is just the usual projection 




(13) 



This inversion formula assumes spherical symmetry. For stacked samples of lenses, 
spherical symmetry follows simply from the statistical isotropy of the Universe. Since we 
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are computing an average mass density of many lenses, it can be written as a correlation 
function, Ap(r) = p & m (r). Equation 12 combined with Equation 11 results in the Equation 
used in Sheldon et al. (2004) to measure the galaxy-mass correlation function: 

a ,n r, , 1 P AE'(i2) + 2AE( J R)AR , 
Ap(r) = fi mPcrii 6m = - / di? 1 7= , -. (14) 

Again, one must truncate this integral since one only has measurements out to some finite 
projected separation R max . We discuss this truncation in Section 5. 



4. Prom Density to Mass Profiles 

Now that we can reconstruct mean 3D density profiles, it would seem straightforward to 
calculate the 3D mass profiles, M(r) = An f^dy y 2 Ap(y). (This mass, the second moment 
of the correlation function, is sometimes referred to as Js(r) in the large-scale structure lit- 
erature, though not usually with cross-correlation functions in mind.) However, the density 
profile cannot simply be integrated directly, since one only has reliable density information 
from shear measurements down to some minimum scale, R m in- One could extrapolate the 
derived density profile to smaller scales, but that would introduce model-dependent assump- 
tions. As we will show, one can do better: the mass inside R m in is in fact constrained by 
the lensing data. This is because Equations 2 and 4 for the tangential shear are non-local. 
To illustrate the point, consider that a point mass creates a shear, 7(7") ~ r~ 2 . However, 
substituting this power-law expression into Equation 5 results in a complete cancellation. 
Thus, adding a point mass does not alter the derived 2D and 3D densities, as should be 
expected. The mass inside R m in is in fact the piece of information that we lost when taking 
the derivative of Equation 4 to get the local Equation 5. Hence, any mass estimator should 
use both the local and non-local equations. 

Throughout this calculation, we will assume that the mean density of the universe has 
been subtracted out, as above. This is just a statement of boundary conditions, which does 
not affect the final 3D profiles, and allows us to ignore the effects of the lensing kernel T, crit . 
We begin with Equation 2 

AE(R) = S(<i?)-E( J R) (15) 

= ^ ( 16 ) 
which defines the mass, m cy i(R), inside a cylinder of radius R oriented along the line of sight: 

m cyl (R) = nR 2 [AZ(R) + E(R)\. (17) 
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Using Equation 13 to replace E(i2), we have 



m cy i(R) = 7T R 2 



POO 

AE(R) +2 drr 

JR 



Ap(r) 



Vr 2 - R 2 



(18) 



Now we want to relate this cylindrical mass, m cy i(R), to the 3D spherical mass M(R). 
It is a slight abuse of notation to write M{R), since R is the 2D radius in the plane of the 
sky, but this avoids the bulkiness of the more precise notation, M(r = R). Again, assuming 
spherical symmetry, the 3D mass can be written as the total mass inside the cylinder minus 
the mass inside the cylinder but outside the sphere. Using cylindrical-polar coordinates 
and spherical radius r = \/x 2 + z 2 , with the z-axis being the axis of projection, we 



x, z 



, A/, 



have 



■ it , ■ 

dx x 2 



7 

J v. 



dz Ap(r) 



M(R) = m cyl (R) - 2tt / 

Jo 

= m cy i(R) - Arc J dr r Ap(r) (r - Vr 2 - R 2 ^ . 



(19) 
(20) 



Substituting Equation 18 for m cy i into 20 gives an expression for M(R) in terms of AS and 
Ap. Since Ap is a functional of AE, this is a mapping: AS i — > M. We will define this 
particular form for M(R) as M out (R), since it only uses data outside of R: 

R 2 



M out (R) = vr/? 2 AE(i?) + 2tt / dr rAp(r) 

'r 



Vr 2 - R 2 



2(r- Vr 2 - i? 2 ) 



(21) 



Since we assumed above that the average density is zero, equation 21 actually gives only the 
mass over-density, not the total mass; to properly correct for this, a term 4/3nR 3 p should be 
added, but this formula is typically used on scales where this extra contribution is negligible. 

We will now derive another formula for the mass profile. As we mentioned, the simplest 
idea is just to integrate the density. We will call this particular form M in (R), which is given 
by 



(22) 
(23) 



rR 

M in (R) = 4tt / dr r 2 Ap(r) 
Jo 

r R 

= M(R min ) + 4tt / dr r 2 Ap(r). 
Replacing the first term with expression 21 evaluated for M out (R min ) gives 

M m {R) = 7cR 2 mm AJ:(R mm ) + 2vr I drr Ap(r) 2r lzJ^ + E(R, R min ), (24) 

** Rmin v ^ min 
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where the last term is 



E(R, R min ) = 2ir 



f 

JR 



dr r Ap(r) 



"min 



'mm 




(25) 



For R ^> Rmin the term E(R, R m i n ) is very small compared to the rest of the terms. Note 
that although we have referred to this expression as M in , the last term does depend on AS 
outside of R. There are thus two seemingly different expressions for the 3D mass profile. In 
fact, one could replace R m in in the above expression with any radius between R m i n and R: all 
of these expressions are equivalent, and they do not add any extra information by combining 
them. They simply add up the information from different radii in a different order, much 
like changing the order of integration on a 2D integral. 

For clusters, one can measure a mean virial mass from these mass profiles quite easily. 
The virial mass is defined as the mass within the virial radius, r V i r , where the virial radius is 
the radius at which the average enclosed density is A vir times either the average or critical 
density of the Universe. Often, A vir = 200 (and critical) is used and this virial radius is 
often denoted r2oo- One must then interpolate between points to find the radius at which 
M(r 2 oo) = 47r/3 r| 00 200 p crit . These non-parametric virial masses should prove very useful 
in comparing data to N-body simulations, which make frequent use of the virial mass. 

As mentioned above, the mass profiles derived from AS actually contain more informa- 
tion than the density profiles. The mass profile at R m i n tells us about the average density at 
R < Rmin, so it can be used to probe these smaller scales. This might be useful in answering 
the question of whether halos have cusps or cores at small radii. It is equally true that 
AS contains all of the information. However projection and non-locality make it somewhat 
harder to use, except through parametric modeling. If one just wants to fit a model to the 
data, one might as well use AS. For visualizing the actual shape of the density profiles 
without a model and obtaining virial mass estimates, these inversions are most useful. 

Finally, we note that there has been some work on deriving aperture mass estimates in 
2D from the shear profile (Schneider 1996; Schneider & Bartelmann 1997). The aperture mass 
m = J d 2 x k(x) w(\x\), where w is an azimuthally symmetric weight function. Schneider 
(1996) suggests choosing a normalized weight function satisfying J d 2 x w{\x\) = 0, i.e., a 
compensated filter, so that the mass estimate is not affected by the mass sheet degeneracy. 
Obviously this would restrict the choice of weight functions so that, for example, the top-hat 
weight function cannot be used. We have implicitly restricted ourselves above to top-hat 
apertures and derived a 3D mass profile from the lensing measurements. We discuss the 
mass sheet degeneracy by showing (in Section 5.2) that, for stacked samples, this unknown 
constant can be predicted fairly well and only creates real uncertainty at the very largest 
scales. 
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5. Practical Application 

With expressions in hand for inverting 3D density and mass profiles from the lensing 
data, AS, we turn to discuss the application of these formulae to real data. The two main 
issues to consider are interpolation to evaluate the various integrals and handling the parts 
of the integrals that extend beyond the scales for which we have measurements. Improper 
handling of the binning and interpolation issues typically results in errors at the few percent 
level for all radii. The issue of extrapolating the integrals to infinity usually only affects the 
recovered profiles at the largest measured scales. 

5.1. Binning and Interpolation 

Typically, one will construct a AS profile by measuring the average weighted lensing 
distortion in a set of radial bins. For example, Sheldon et al. (2004) measured AS(i?) in 
18 logarithmically-spaced radial bins from 0.02 hr l Mpc to 11 hr x Mpc, giving a bin width 
Ri + i/Ri = 1.41. One would like to assign one effective radius to each bin. The radius 
assigned should be that radius at which we believe the true AS(i?) is equal to our measured 
average. The mean value theorem of integral calculus guarantees that this radius will be 
inside the bin if AS(i?) is continuous, yet it is impossible to find this correct value precisely 
without first knowing the shape of the function within the bin. In practice, however, it is 
good enough to approximate the function as a local power-law, for which the correct radius 
can be calculated exactly. The power-law exponent should roughly match the slope of the 
data though the result is not very sensitive to the exact value. This should be the first step 
carried out before applying the density and mass formulas. For the binning and data of 
Sheldon et al. this procedure produces only a ~ 1% shift in effective radius from simply 
using the midpoint, mean, or geometric mean of the bin boundaries Ri and Ri+\. For larger 
radial bins, say, as large as R i+ i/Ri = 1.8, this interpolation can lead to a 5% shift in the 
assigned radius. 

After selecting the proper effective radii, Ri, corresponding to each AS(i?j) bin, one can 
apply the inversion formulae of §4. This requires a scheme for interpolating AS between the 
bins and for interpolating Ap between the points to infer the mass. Saunders et al. (1992) 
applied Von Zeipel's inversion formula, essentially our Equation 12, to the projected galaxy 
autocorrelation function, w p (R), to obtain the 3D autocorrelation function, £(r). They used 
a piecewise-linear model for w p (R) and hence a piecewise-constant model for w'(R), which 
enters into the Von Zeipel integral. This zero-th order interpolation allows calculation of 
the integral over each bin analytically, which results in simply a weighted sum over the data 
points (details can be found in Saunders et al. (1992)). However, the data for both w p and 
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AE look roughly like power laws, AS ~ R~ - 8 , so this simple interpolation scheme biases 
the inversion low by about 10% for Sheldon et al.-sized bins. If one takes this scheme to 
first order by representing the derivative, AE'(i2), as piece-wise linear rather than constant 
over a bin, this biases the result by about the same amount, but high instead of low. One 
can in fact do n-point Lagrangian interpolation to any order and still evaluate the integrals 
analytically over each bin. Such methods are completely linear in the data, which is another 
advantage. In principle, higher order interpolation will reduce this bias, but only at the 
expense of correlating many neighboring points. We have found in practice that the best 
method is to use a power-law interpolation, equivalent to linear interpolation of the log- 
log values. Any two positive data points define a unique power-law between them. The 
integrals for each bin need to be evaluated numerically, so this is slower computationally, 
and care needs to be taken to handle negative data points which may be present for noisy 
data. We therefore do not take logarithms directly but rather model the function as a power 
law plus constant, where the constant is only used for negative values. Data for which 
there are negative values are generally very noisy; for that case, interpolation is never the 
dominant source of errors. Tests of the inversion show that for power-law density profiles, 
this interpolation scheme is exact (unsurprisingly); for NFW profiles and halo- model type 
profiles, the power-law interpolation is good to a fraction of a percent for Sheldon et al.-sized 
bins (i.e. R i+ \/ Ri = 1.41) and very robust for noisy data. 



5.2. Endpoint Corrections 

The inversion integrals in Sections 3 and 4 should in principle be evaluated from the 
minimum radial bin, R m i n , to oo. In practice, one only has shear data covering a range 
from R m in to R m ax- We refer to the parts of the integrals from R max to oo as endpoint 
corrections; they must be added to the parts of the integrals that we can actually perform 
by interpolating the data. Consider Equation 12 as an example; we can rewrite it as 

Ap(r) = I/ dR ^B= + 1 iR ^£L^T D+ T c , (26) 

J r VR 2 -r 2 nJ Rmax y/R 2 -r 2 

where T D and T c stand for Data and Correction. To estimate the magnitude of each term, 
consider an isothermal sphere, AS ~ R^ 1 ', in this case, the integrals can be computed 
analytically, leading to 

Td + Tc \j ( Rmax ) ^ ^ 

For the outermost data point, r = R max , the correction term gives the whole result, and 
one cannot learn anything about the density at this point. However, for smaller scales, 



-14- 



r <C Rmax, the ratio T C /(T D +T C ) ^ (r / R max ) 2 /2. For example, for Sheldon et al.-size bins, 
this implies that the correction is only about 14% at the value of Ri = R±q, two radial bins 
inward of R m ax — Ris- The correction terms at the next bins inward (proceeding from larger 
to smaller scale) are 7% , 3%, and 2%, while the remaining inner 12 bins are affected by less 
than 1%. In fact, it is only the error in estimating this correction term that concerns us, so 
we expect that even at the second bin inward of R ma x, one could estimate Ap to an accuracy 
of a few percent. Around the virial radius of clusters, R ~ 1 h^ 1 Mpc, the correction term 
itself is only half a percent of the total for R max = 10 hr x Mpc. 

Beyond 10/i -1 Mpc, the halo density profile is expected to drop off more steeply than 
isothermal, so the correction factor in this case should be less than the estimate in 27. At 
these large scales, it is probably a good assumption to use the shape of the linear mass 
correlation function (i.e. the Fourier transform of the linear power spectrum), which is well 
understood theoretically for a given cosmology, as the model with which to compute the 
endpoint correction. This is the expected large scale behavior from the halo model. This 
requirement of providing endpoint corrections is the one part of this method which requires a 
model. However, this model dependence only affects the results at the largest scales for any 
reasonable extrapolation, guided by theoretical expectation and constrained by the large- 
scale galaxy autocorrelation function which is now measured to 175 br x Mpc (Eisenstein 
et al. 2005). It is left up to the experimenter to decide which points to omit and how to take 
the error in the correction into account. 

In this context, for lensing data that extends to scales of order a few Mpc, the endpoint 
correction can be considered a manifestation of the mass sheet degeneracy. In the current 
paradigm of structure formation, the Universe is not filled with uniform-density mass sheets, 
but it does contain large-scale density perturbations that are correlated with dark matter 
halos. For a cluster positioned near the peak of a large-scale over-density with scale 3> R max , 
the latter can be thought of to first approximation as a mass sheet for the cluster shear 
measurement. Using the theoretical mass correlation function on large scales to estimate the 
endpoint correction, as described above, therefore provides an estimate of the effective mass 
sheet degeneracy for this method. As we have argued above and show below in an N-body 
simulation, the effect of this degeneracy on the estimate of cluster virial masses is below 1%. 



5.3. Propagation of Errors 

In performing these inversions, it is important to correctly propagate the errors from 
AE(i?) to Ap(r) and M(r). In particular, the covariance matrix between radial bins Ri needs 
to be calculated. We will assume that we have already calculated the covariance matrix for 
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AEj which we denote by C A . Usually this is just diagonal from shot noise (i.e., lensing 
shape noise) but at larger scales there can be off-diagonal terms from sample variance and 
systematics such as imperfect PSF anisotropy correction. These off-diagonal terms can be 
estimated with jackknife resampling (Sheldon et al. 2004). We will denote the covariance 
matrix for Api and Mj as C?- and . 

As we have mentioned, with the simplest polynomial interpolations, the mapping AE(i?) i — > 
Ap(r) is linear, so there exists a matrix A such that Api = A^-AE^; likewise, since M is 
a linear combination of AE and Ap, and Ap is computed linearly from AS, there exists a 
matrix B such that M{ = BijAHj. For the polynomial interpolation schemes, these matrix 
elements can be computed analytically, without requiring numerical integration. For the 
power-law interpolation that we recommend using, the mapping is not exactly linear, but it 
is well approximated by its first order Taylor expansion, 

dAp^ jm_ 

13 dAZj ' 13 i)A\, 

Although, these derivatives can be calculated analytically in principle, in practice we estimate 
them numerically with a finite difference method. 

Alternatively, one could use the linear interpolations just for the covariance matrices, 
where the matrix elements can be computed analytically. This would be the fastest method 
computationally. Once the linearization has been computed (by any method), the covariance 
matrices can be propagated simply with matrix multiplication, C p = A C A A T and C M = 
B C A B T . Generally, we find that the covariance matrix C p is mostly diagonal if C is 
diagonal. The covariance matrix C M , however, has large off-diagonal terms since M(r) is a 
cumulative statistic and therefore neighboring bins are correlated. 

6. Tests of the Method: N-body Simulations 

As a proof of principle, we have performed tests of these inversion methods on an N- 
body CDM simulation. In the simulation, we can measure the 3D density and mass profiles 
directly and check that the inversion of projected quantities correctly recovers the true values. 

The simulation we use has 512 3 particles in a periodic cube of length 300 h^ 1 Mpc. The 
simulation is evolved from z = 60 to z = using a TreePM code (White 2002, 2003) 1 ; 
we use only the z = output. The cosmological parameters used are Qm = 0.3, = 0.7, 



The simulation is publicly available at http:/ /mwhite. berkeley.edu/Siml/ 



(28) 
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h = 0.7, n — 1, VL b h 2 = 0.02, and <t 8 = 1. The simulation has an effective Plummer force- 
softening scale of 20 h' 1 kpc which is fixed in comoving coordinates. The mass of each dark 
matter particle is 1.7 x lO lo /i _1 M . Dark matter halos are identified using a Friends-of- 
Friends (FoF) algorithm (Davis et al. 1985) with a linking length of 0.2 in units of the mean 
inter-particle separation. Specific details of the simulation such as resolution, cosmology, 
and halo finding are not crucially important, since we are only interested in whether the 
inversion methods recover the 3D quantities. The simulation box was chosen to have high 
enough resolution to resolve the inner regions of clusters measurable by SDSS and is large 
enough to have relatively low cosmic variance; nevertheless, it is quite a bit smaller than the 
size of the final SDSS cluster sample currently in preparation. 

We select all halos of mass M vir > lO 14 ^" 1 Mpc, measure p(r) and M(r), and average 
these 3D quantities for all such massive halos. For p(r), we correct for the effects of bin- 
ning in a similar way to the methods described in Section 5.1. Next, we project the box 
separately along each of its three axes (x,y,z). For each of these projections, we measure the 
average E(i?) and AE(i?) from Equation 2. We do not perform ray tracing to determine the 
real shear, so we are only testing the inversion method, assuming that one can measure AE 
through measurement of background galaxy shear (see Sheldon, et al.). In particular, we are 
not testing the shear non-linearities, the details of galaxy shape measurement, photometric 
redhshift estimation and calibration, or any kind of cluster finding. These details are impor- 
tant for any analysis of real data, but they are beyond the scope of this paper (and will be 
discussed further in future papers in this series). 

By studying the differences in these projected 2D quantities for the three different 
projections of the simulation box, we can crudely measure a kind of sample variance which 
we call anisotropic sample variance (ASV). If the halo- mass correlation function were exactly 
isotropic, all three projections would agree precisely. Since this is only true in the infinite 
volume limit, there will be differences between the three projections of the simulation as well 
as small departures from isotropy in any real measurement. On small scales, these differences 
arise from the asphericity of halos: if halo orientations are not strongly correlated, the 
residual ellipticity will result in a small random error that will decrease as ~ 1/ V Nhaios ■ More 
explicitly, we can define the halo ellipticity along the ^-direction as e z = (3a 2 — a 2 ) /a 2 , where 
a 2 is the second moment of the cluster density distribution along the projected direction and 
a 2 = a 2 + a 2 + a 2 defines the characteristic cluster size. 

The ensemble average of e z is zero by symmetry, but the rms for a finite sample will 
be zero within ±^/< e 2 > /Nhaios- One can show that this rms halo ellipticity induces a 
multiplicative correction factor in the projected density and all derived quantities (p and 
M) of 1 ± a/< e 2 > ]N r halos] this factor will vary randomly from realization to realization (or 
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from projection to projection within a single N-body realization). The rms halo ellipticity 
a/ < e 2 > appears to be about 0.5 for the massive halos in the N-body simulation described 
above. Since large survey samples will have very many clusters, this correction will usually 
be a very small effect. For example, in the sample used in this study, there are 1226 halos 
above the mass threshold in the simulation, resulting in an expected multiplicative variance 
of ±1.4% in the inferred profiles. This estimate agrees very well with the observed scatter 
between the three different projections on small and intermediate scales, as shown below in 
Figures 2 and 3. The SDSS cluster samples are typically larger than this, so the effect will 
be correspondingly reduced. Less massive halos (these are M = lO 14 /i _1 M ), being more 
common, will have a larger N ha [ os and thus this will be reduced further. 

On larger scales, where the halo-mass correlation function is dominated by contributions 
from mass elements in halos different from the lens halo (i.e., where, in halo model terms, the 
two-halo term dominates (Seljak 2000; Mandelbaum et al. 2004)), the asymmetric sample 
variance arises predominantly from the shapes of larger-scale structures such as filaments 
and superclusters. Since there are fewer structures contributing to this part of the ASV 
(fewer filaments than halos), and since these large structures are more asymmetric than the 
halos, the variance between projections is expected to be larger on large scales than small 
scales, in qualitative agreement with Figure 2 below. This large-scale asymmetric sample 
variance is expected to scale as 1/ VV where V is the volume of the survey. When the shot 
noise from galaxy shapes (i.e., the shape noise) is sufficiently small, the large-scale errors will 
be dominated by the ASV. Judging by the volume of our simulation and the scatter that we 
see at scales of a few Mpc/h, we expect this large-scale multiplicative error to scale roughly 
as 2% x (V/Gpc 3 )~ - 5 , although it may vary with scale R. 

Figure 1 shows the AS results for each of the three projections of the simulation (the 
three different colors correspond to the three different projections). We note that these curves 
turn over below about 20/i _1 kpc, an effect almost certainly due to the finite resolution of the 
simulation: this scale is about 10 times the Plummer force-softening scale. The turn-over 
in AS typically occurs at a few times the 3D correlation function core radius and the 3D 
correlation function core radius typically occurs at a few times the Plummer force softening 
scale. By correlation function core radius we have in mind something like a power-law with 
a softening, r — > ^/r 2 + r 2 ore . Smaller-volume, higher-resolution simulations do not show 
this turn-over in AS. For a pure NFW profile, AS flattens at small scale but does not 
decline. As noted above, despite this resolution issue, we should be able to use these small 
scales to test our inversion methods since the 3D values are equally affected by resolution. 
In fact, this gives us an opportunity to test the robustness of the methods by including a 
pathological regime instead of a perfect power-law or NFW profile on all scales. The large 
errors (estimated by a jackknife technique) at small scales indicate that there are not many 
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particles in these weak cusps. The larger deviations at the largest scales indicate the larger 
ASV, as expected from the arguments above. 

We invert each of these AE profiles to obtain the estimated density Ap(r) and M(r) 
using the formulae and algorithms described in previous sections. The endpoint corrections 
are carried out assuming power-law extrapolations with logarithmic slopes that are not too 
different from the expected continuation of the linear theory correlation function. The last 
two bins from the inversions are not plotted, since they are highly affected by the endpoint 
corrections. 

Figure 2 shows results of the inversion for the density profile p(r). As in Fig. 1, the 
red, green, and blue curves denote the inversions from the three different projections. The 
magenta line shows the true 3D density as measured in spherical shells. On this log-log plot 
one can hardly see any difference between the true and inverted profiles. The bottom panel 
of Figure 2 shows the ratio of each of these inverted density profiles to the true 3D density 
with a linear axis. On scales r < 1 h^ 1 Mpc, these ratios are all consistent with unity to 
within about 5%. On larger scales, the ASV becomes larger, and the three curves differ 
randomly from unity by about 20%. There do not appear to be any systematic biases or 
significant trends. Changing the endpoint correction moves these curves around by a few 
percent at the largest scales. 

Figure 3 shows results of the inversions for the mass profile, M(r). Again, the top panel 
shows all three inversions with the true M(r) over-plotted in magenta; the inversions and the 
true profile are indistinguishable on this plot. The lower panel shows the ratios of inverted 
to true mass. Except for the smallest scales, where the scatter between the three inversions 
is a little larger, the mass ratios are consistent with unity to about 5%. Around the virial 
radius, r ~ lh^ 1 Mpc, which is the optimal location in terms of measurement signal-to- 
noise, the inverted masses are correct to about 2%. There appears to exist a slight tilt to 
these ratio curves, but the points at different radii are strongly correlated (M(r) being a 
cumulative quantity), so this trend is not significant. On the whole, we consider these tests 
a confirmation that inverting lensing measurements can correctly recover the density and 
mass profiles of galaxy and cluster halos. 

7. Effects of Halo Asphericity and Large-scale Structure along the Line of 

Sight 



In the case of weak lensing mass reconstructions for individual clusters, one must con- 
sider the systematic effects of both cluster asphericity and contamination from structure 
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Fig. 1. — The mean tangential shear A£(i?) for halos with virial masses greater than 
lO 14 /i _1 M for three orthogonal projections of a ACDM N-body simulation. The turn- 
over below R ~ 20 h~ x kpc is due to the finite resolution of the simulation. The scatter 
between the curves at the smallest scale is likely due to shot noise since there are very few 
dark matter particles at these scales due to the resolution. At large radius, we attribute the 
differences to anisotropic sample variance (ASV) caused by large-scale structure along the 
line of sight. 
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Fig. 2. — Top panel: Inverted 3D density profiles for the three projections of the simulation. 
The red , green, and blue curves indicate the inversions, while the magenta curve denotes 
the true density profile measured by averaging the dark matter particle density in spherical 
shells. The inverted profiles are only distinguishable from each other and from the true 
profile at large scales. Bottom panel: ratio of the inverted density to the true 3D density 
for the three simulation projections as a function of radius; curves correspond to the same 
projections as in the top panel. For each projection, the inverted and true density profiles 
agree to about 5% or better within the virial radius (r 2 oo = 0.585 hr 1 Mpc). On larger 
scales, there is significantly more scatter due to ASV in this relatively small simulation box. 
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Fig. 3. — Top panel: mass profiles inverted from the tangential shear for the three projections 
of the N-body simulation (red, green, and blue curves). The magenta line shows the true 
average mass profile measured in 3D for these halos. The differences are not evident on 
this scale. Bottom panel: ratios of the three inverted mass profiles to the true 3D profile. 
The inverted masses are all correct to within 5% out to the largest scales and to within 2% 
near the virial radius. Because the points are correlated, the slight trend with radius that 
is apparent is not significant. The deviations at the smallest scale are mostly likely due to 
shot noise. 
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along the line of sight, such as filaments or clusters at other distances seen in projection. 
Several studies (Cen 1997; Metzler et al. 1999, 2001; White et al. 2002) have concluded 
that large-scale structure (LSS) along the line of sight introduces significant scatter and bias 
in the resulting cluster mass estimates. Later studies (Hoekstra 2001, 2003; Clowe et al. 
2004; de Putter & White 2004) agreed that LSS introduces scatter in the mass estimate, 
though Clowe et al. (2004) argued that the bias found in the previous studies was due to an 
over-simplification of the mass estimation. Dodelson (2004) suggested ways that the extra 
noise from LSS can be reduced, though not eliminated, through a particular type of filtering 
when making mass maps. Clowe et al. (2004), on the other hand, argued that cluster halo 
asphericity is a larger problem for determining cluster masses than LSS along the line of 
sight: if the major (minor) axis of the halo is aligned along the line of sight, the virial mass 
will be overestimated (underestimated). In general one does not know the orientation of the 
halo in 3D, so this creates an uncertainty in the inferred mass at the ~ 30% level. 

As demonstrated in Section 6 and Figure 3, neither of these effects compromises the 
accuracy of our statistical mass inversion method, provided the cluster sample is sufficiently 
large. By stacking many clusters to form a statistical sample, the effects of large-scale 
structure uncorrelated with the lensing cluster must average out by statistical isotropy. On 
small scales, the mean orientation of aspherical halos for a large sample of lenses must be 
very close to random, provided that the method of selecting cluster lenses is not biased 
toward including those oriented, say, along the line of sight. These selection effects can be 
checked in N-body simulations that incorporate some prescription for assigning luminous 
galaxies to dark matter halos; as such, they are beyond the scope of this paper, but we 
plan to address them in a future publication. (Note that cluster samples selected by their 
shear signal may suffer from such a bias (White et al. 2002).) On intermediate scales, the 
orientation of nearby filamentary structure correlated with the cluster lens (as well as the 
angular distribution of nearby correlated halos) must also be random for a large statistical 
sample. Since such structure is correlated with the lens, it does contribute to the tangential 
shear on large scales; indeed, as we have argued above, it dominates the signal from clusters 
on scales above a few /i~ 1 Mpc and simply represents the large-scale cluster-mass correlation 
function. Again, as shown above, the effects of correlated structure along the line of sight 
appear to be negligible out to the virial radius of massive clusters. 

The only residual problem for the statistical inversion method is anisotropic sample 
variance at large scales (discussed above in Section 6); this finite- volume effect creates a 
small scatter in the inferred mass but does not create a bias. 



-23- 



8. Discussion 

We have developed a non-parametric method for inverting cross-correlation lensing mea- 
surements to obtain average density and mass profiles of halos. We have demonstrated the 
method on an N-body simulation and shown that it successfully recovers the 3D profiles. 
We argued that asphericity of halos and large-scale structure along the line of sight do not 
introduce bias or substantial uncertainty in the inferred mass estimates. 

This method should find several useful applications in surveys. Applied to galaxy- 
galaxy lensing measurements, this method can be used to measure Ap(r) = Q m p crit £ sm (r) for 
samples of lens galaxies (e.g., Sheldon, et al. 2004). On small scales, this provides information 
about galaxy dark matter halos and should thereby constrain models of galaxy formation and 
evolution, including the details of hierarchical structure formation and merging. Statistical 
mass profiles around galaxies can provide virial mass measurements and, combined with the 
average light profile around these galaxies, determine mass-to-light ratios as a function of 
scale. 

On larger scales, galaxy-galaxy lensing inversions should tell us more about cosmology. 
The autocorrelation function of galaxies allows us to measure ^ gg {r) = b 2 (r)^ mm (r), where 
b(r) is the scale-dependent bias and £ mm (r) is the mass autocorrelation function. Lensing 
allows us to measure Q m ^ gm {r) = Q m b(r) r x (r) £ mm (r), where the cross-bias, r x (r), is 
sometimes referred to as stochastic bias (a term we do not recommend using in this context, 
since it is not bounded by ±1). The presence of these two bias functions, b(r) and r x (r), 
allows for the most general model of the relative clustering of galaxies and mass. (See also 
Neyrinck et al. (2005).) Simulations indicate that r x (r) is consistent with unity on scales 
larger than ~ 2h~ l Mpc for any range of halo masses (Tasitsiomi et al. 2004). In this case, 
measurement of large-scale galaxy-galaxy lensing and of the galaxy autocorrelation function 
determine the mass correlation function, f2^£ mm = (fl m £, g m) 2 / (, gg and so fixing the shape 
of the linear correlation function determines the amplitude or normaliztion Q m cr s . It also 
determines the bias, b(r)/Q m = ^ gg /^i m ^ g m- The assumption that bias is scale independent on 
large scales is crucial to extracting cosmological information from the galaxy power spectrum; 
lensing can provide a way to test this important assumption. Other cosmological probes, 
such as cluster counts and cosmic shear, measure a different parameter combination, fl^as, 
so combining these probes with cross-correlation lensing helps break this classic degeneracy. 
An advantage of this new method is that it is robust; the only assumptions used are that 
general relativity is correct and that r x (r) ~ 1 on large scales. 

Applied to galaxy clusters, this inversion method should prove useful for both astro- 
physics and cosmology. On small scales, cross-correlation lensing probes the mean density 
profiles of clusters. Cosmological dark matter simulations indicate that clusters have univer- 
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sal dark matter halos (though baryonic physics presumably needs to be taken into account 
to understand the inner density structure in detail). Inverted density profiles can test this 
assumption quite directly. As we have demonstrated, cluster mass profiles can be deter- 
mined, leading to virial mass estimates independent of a model for the density profiles. 
Virial masses allow a direct comparison between simulated and real clusters. In addition, 
in a cluster survey we can measure the cluster abundance, n(M,z), binned by any observ- 
able proxy for mass, and use this inversion method to calibrate the virial mass-observable 
relation as a function of redshift. This should allow cluster surveys to more precisely probe 
cosmology, including the dark energy. We still need to gather information about the scatter 
in the mass-observable relation; this can be constrained by simulations, by self-calibration 
(Lima & Hu 2005), and internally by dividing the clusters into subsamples and separately 
estimating the mass-observable relation. 

Before applying this method with confidence to real cluster data, the possible selection 
biases involved in a given cluster selection algorithm need to be explored. If the selection 
algorithm preferentially finds clusters aligned along the line of sight, the assumption of 
statistical isotropy will be violated and the mass inversions may be biased. Similarly, if the 
selection algorithm assigns cluster centroids that are displaced from the true halo centers, the 
resulting density profiles will be convolved with distribution of centroid errors. Additional 
complications for the cluster counting technique arise if clusters, however selected, are not 
isomorphic to massive dark matter halos. These issues can be tested with simulations. 

Applying these lens inversion methods to clusters on large scales provides new cosmolog- 
ical information beyond the cluster counting technique. As noted for galaxies above, we can 
use the cluster-mass and cluster-cluster correlations to measure fl m (Js- A useful cross-check 
is provided by comparing the results for different cluster samples binned by some observable, 
e.g., richness. This method should be even more robust for clusters than for galaxies, since 
the halo bias as a function of mass is predictable from theory and can be studies in dissi- 
pationless N-body simulations (Seljak & Warren 2004; Tasitsiomi et al. 2004). One can use 
this information in two ways. Conservatively, one could measure the average cluster masses, 
M, and b(M)/Q m ; from these bias measurements one could further constrain cosmology 
through the theoretical bias predictions and the direct measurement of fi m cr 8 . Alternatively, 
one could use the lensing data only on small scales to determine masses and assume that 
the theoretical bias predictions are correct; in this approach, one would use the mass and 
an assumed cosmology to predict the bias. This would allow one to 'de-bias' the cluster 
autocorrelation function and obtain a constraint on a 8 that depends on other cosmological 
parameters as well. Since the scaling of bias with mass is somewhat weak (b(M) ~ M a2 in 
the vicinity of M = 1O 13 /i _1 M ), and the cluster autocorrelation function can be measured 
with high signal-to-noise, this measurement of Q m os could be quite precise. This method 
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has been applied to galaxies by Seljak et al. (2005), utilizing the halo model. We believe that 
clusters are a more natural choice for this method in that they do not need an additional 
halo occupancy presciption. 

The approach of 'de-biasing' the correlation function also makes possible the study of 
the growth of structure. In principle, the inversion method can be applied to data in redshift 
bins, allowing one to measure the large-scale mass power spectrum amplitude as a function 
of redshift. This is particularly useful, because the linear density perturbation growth rate 
is sensitive to dark energy. Moreover, since the evolution of the halo mass function is also 
sensitive to the growth of structure, combining cluster counts and the cluster-mass correlation 
function should offer a powerful probe of cosmology. 

We are currently applying these statistical lensing inversion techniques to optically se- 
lected clusters in the Sloan Digital Sky Survey. In the future, on-going and planned wide-area 
surveys, including RCS II, CFHTLS, the VST and VISTA surveys, the Dark Energy Survey, 
PanSTARRS, LSST, and SNAP/JDEM will allow for higher signal-to-noise lensing mea- 
surements to be made and will eventually probe much larger cosmological volumes. The 
cosmological constraints they supply from cross-correlation lensing will complement those 
from other probes. 
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